!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief provides a resp fit for gas phase systems
!> \par History
!>      created
!>      Dorothea Golze [06.2012] (1) extension to periodic systems
!>                               (2) re-structured the code
!> \author Joost VandeVondele (02.2007)
! **************************************************************************************************
MODULE qs_resp
   USE atomic_charges,                  ONLY: print_atomic_charges
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind
   USE bibliography,                    ONLY: Campana2009,&
                                              Golze2015,&
                                              Rappe1992,&
                                              cite_reference
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell,&
                                              pbc,&
                                              use_perd_none,&
                                              use_perd_xyz
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_generate_filename,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
   USE cp_units,                        ONLY: cp_unit_from_cp2k,&
                                              cp_unit_to_cp2k
   USE input_constants,                 ONLY: do_resp_minus_x_dir,&
                                              do_resp_minus_y_dir,&
                                              do_resp_minus_z_dir,&
                                              do_resp_x_dir,&
                                              do_resp_y_dir,&
                                              do_resp_z_dir,&
                                              use_cambridge_vdw_radii,&
                                              use_uff_vdw_radii
   USE input_section_types,             ONLY: section_get_ivals,&
                                              section_get_lval,&
                                              section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kahan_sum,                       ONLY: accurate_sum
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE machine,                         ONLY: m_flush
   USE mathconstants,                   ONLY: pi
   USE memory_utilities,                ONLY: reallocate
   USE message_passing,                 ONLY: mp_para_env_type,&
                                              mp_request_type
   USE particle_list_types,             ONLY: particle_list_type
   USE particle_types,                  ONLY: particle_type
   USE periodic_table,                  ONLY: get_ptable_info
   USE pw_env_types,                    ONLY: pw_env_get,&
                                              pw_env_type
   USE pw_methods,                      ONLY: pw_copy,&
                                              pw_scale,&
                                              pw_transfer,&
                                              pw_zero
   USE pw_poisson_methods,              ONLY: pw_poisson_solve
   USE pw_poisson_types,                ONLY: pw_poisson_type
   USE pw_pool_types,                   ONLY: pw_pool_type
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r3d_rs_type
   USE qs_collocate_density,            ONLY: calculate_rho_resp_all,&
                                              calculate_rho_resp_single
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type,&
                                              set_qs_env
   USE qs_kind_types,                   ONLY: qs_kind_type
   USE qs_subsys_types,                 ONLY: qs_subsys_get,&
                                              qs_subsys_type
   USE uff_vdw_radii_table,             ONLY: get_uff_vdw_radius
#include "./base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

! *** Global parameters ***

   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_resp'

   PUBLIC :: resp_fit

   TYPE resp_type
      LOGICAL                                :: equal_charges = .FALSE., itc = .FALSE., &
                                                molecular_sys = .FALSE., rheavies = .FALSE., &
                                                use_repeat_method = .FALSE.
      INTEGER                                :: nres = -1, ncons = -1, &
                                                nrest_sec = -1, ncons_sec = -1, &
                                                npoints = -1, stride(3) = -1, my_fit = -1, &
                                                npoints_proc = -1, &
                                                auto_vdw_radii_table = -1
      INTEGER, DIMENSION(:), POINTER         :: atom_surf_list => NULL()
      INTEGER, DIMENSION(:, :), POINTER       :: fitpoints => NULL()
      REAL(KIND=dp)                          :: rheavies_strength = -1.0_dp, &
                                                length = -1.0_dp, eta = -1.0_dp, &
                                                sum_vhartree = -1.0_dp, offset = -1.0_dp
      REAL(KIND=dp), DIMENSION(3)            :: box_hi = -1.0_dp, box_low = -1.0_dp
      REAL(KIND=dp), DIMENSION(:), POINTER   :: rmin_kind => NULL(), &
                                                rmax_kind => NULL()
      REAL(KIND=dp), DIMENSION(:), POINTER   :: range_surf => NULL()
      REAL(KIND=dp), DIMENSION(:), POINTER   :: rhs => NULL()
      REAL(KIND=dp), DIMENSION(:), POINTER   :: sum_vpot => NULL()
      REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix => NULL()
   END TYPE resp_type

   TYPE resp_p_type
      TYPE(resp_type), POINTER              ::  p_resp => NULL()
   END TYPE resp_p_type

CONTAINS

! **************************************************************************************************
!> \brief performs resp fit and generates RESP charges
!> \param qs_env the qs environment
! **************************************************************************************************
   SUBROUTINE resp_fit(qs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env

      CHARACTER(len=*), PARAMETER                        :: routineN = 'resp_fit'

      INTEGER                                            :: handle, info, my_per, natom, nvar, &
                                                            output_unit
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ipiv
      LOGICAL                                            :: has_resp
      REAL(KIND=dp), DIMENSION(:), POINTER               :: rhs_to_save
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(section_vals_type), POINTER                   :: cons_section, input, poisson_section, &
                                                            resp_section, rest_section

      CALL timeset(routineN, handle)

      NULLIFY (logger, atomic_kind_set, cell, subsys, particles, particle_set, input, &
               resp_section, cons_section, rest_section, poisson_section, resp_env, rep_sys)

      CPASSERT(ASSOCIATED(qs_env))

      CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, input=input, &
                      subsys=subsys, particle_set=particle_set, cell=cell)
      resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
      CALL section_vals_get(resp_section, explicit=has_resp)

      IF (has_resp) THEN
         logger => cp_get_default_logger()
         poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
         CALL section_vals_val_get(poisson_section, "PERIODIC", i_val=my_per)
         CALL create_resp_type(resp_env, rep_sys)
         !initialize the RESP fitting, get all the keywords
         CALL init_resp(resp_env, rep_sys, subsys, atomic_kind_set, &
                        cell, resp_section, cons_section, rest_section)

         !print info
         CALL print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)

         CALL qs_subsys_get(subsys, particles=particles)
         natom = particles%n_els
         nvar = natom + resp_env%ncons

         CALL resp_allocate(resp_env, natom, nvar)
         ALLOCATE (ipiv(nvar))
         ipiv = 0

         ! calculate the matrix and the vector rhs
         SELECT CASE (my_per)
         CASE (use_perd_none)
            CALL calc_resp_matrix_nonper(qs_env, resp_env, atomic_kind_set, particles, cell, &
                                         resp_env%matrix, resp_env%rhs, natom)
         CASE (use_perd_xyz)
            CALL cite_reference(Golze2015)
            IF (resp_env%use_repeat_method) CALL cite_reference(Campana2009)
            CALL calc_resp_matrix_periodic(qs_env, resp_env, rep_sys, particles, cell, natom)
         CASE DEFAULT
            CALL cp_abort(__LOCATION__, &
                          "RESP charges only implemented for nonperiodic systems"// &
                          " or XYZ periodicity!")
         END SELECT

         output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
                                            extension=".resp")
         IF (output_unit > 0) THEN
            WRITE (output_unit, '(T3,A,T69,I12)') "Number of fitting points "// &
               "found: ", resp_env%npoints
            WRITE (output_unit, '()')
         END IF

         !adding restraints and constraints
         CALL add_restraints_and_constraints(qs_env, resp_env, rest_section, &
                                             subsys, natom, cons_section, particle_set)

         !solve system for the values of the charges and the lagrangian multipliers
         CALL DGETRF(nvar, nvar, resp_env%matrix, nvar, ipiv, info)
         CPASSERT(info == 0)

         CALL DGETRS('N', nvar, 1, resp_env%matrix, nvar, ipiv, resp_env%rhs, nvar, info)
         CPASSERT(info == 0)

         IF (resp_env%use_repeat_method) resp_env%offset = resp_env%rhs(natom + 1)
         CALL print_resp_charges(qs_env, resp_env, output_unit, natom)
         CALL print_fitting_points(qs_env, resp_env)
         CALL print_pot_from_resp_charges(qs_env, resp_env, particles, natom, output_unit)

         ! In case of density functional embedding we need to save  the charges to qs_env
         NULLIFY (dft_control)
         CALL get_qs_env(qs_env, dft_control=dft_control)
         IF (dft_control%qs_control%ref_embed_subsys) THEN
            ALLOCATE (rhs_to_save(SIZE(resp_env%rhs)))
            rhs_to_save = resp_env%rhs
            CALL set_qs_env(qs_env, rhs=rhs_to_save)
         END IF

         DEALLOCATE (ipiv)
         CALL resp_dealloc(resp_env, rep_sys)
         CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
                                           "PRINT%PROGRAM_RUN_INFO")

      END IF

      CALL timestop(handle)

   END SUBROUTINE resp_fit

! **************************************************************************************************
!> \brief creates the resp_type structure
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
! **************************************************************************************************
   SUBROUTINE create_resp_type(resp_env, rep_sys)
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys

      IF (ASSOCIATED(resp_env)) CALL resp_dealloc(resp_env, rep_sys)
      ALLOCATE (resp_env)

      NULLIFY (resp_env%matrix, &
               resp_env%fitpoints, &
               resp_env%rmin_kind, &
               resp_env%rmax_kind, &
               resp_env%rhs, &
               resp_env%sum_vpot)

      resp_env%equal_charges = .FALSE.
      resp_env%itc = .FALSE.
      resp_env%molecular_sys = .FALSE.
      resp_env%rheavies = .FALSE.
      resp_env%use_repeat_method = .FALSE.

      resp_env%box_hi = 0.0_dp
      resp_env%box_low = 0.0_dp

      resp_env%ncons = 0
      resp_env%ncons_sec = 0
      resp_env%nres = 0
      resp_env%nrest_sec = 0
      resp_env%npoints = 0
      resp_env%npoints_proc = 0
      resp_env%auto_vdw_radii_table = use_cambridge_vdw_radii

   END SUBROUTINE create_resp_type

! **************************************************************************************************
!> \brief allocates the resp
!> \param resp_env the resp environment
!> \param natom ...
!> \param nvar ...
! **************************************************************************************************
   SUBROUTINE resp_allocate(resp_env, natom, nvar)
      TYPE(resp_type), POINTER                           :: resp_env
      INTEGER, INTENT(IN)                                :: natom, nvar

      IF (.NOT. ASSOCIATED(resp_env%matrix)) THEN
         ALLOCATE (resp_env%matrix(nvar, nvar))
      END IF
      IF (.NOT. ASSOCIATED(resp_env%rhs)) THEN
         ALLOCATE (resp_env%rhs(nvar))
      END IF
      IF (.NOT. ASSOCIATED(resp_env%sum_vpot)) THEN
         ALLOCATE (resp_env%sum_vpot(natom))
      END IF
      resp_env%matrix = 0.0_dp
      resp_env%rhs = 0.0_dp
      resp_env%sum_vpot = 0.0_dp

   END SUBROUTINE resp_allocate

! **************************************************************************************************
!> \brief deallocates the resp_type structure
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
! **************************************************************************************************
   SUBROUTINE resp_dealloc(resp_env, rep_sys)
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys

      INTEGER                                            :: i

      IF (ASSOCIATED(resp_env)) THEN
         IF (ASSOCIATED(resp_env%matrix)) THEN
            DEALLOCATE (resp_env%matrix)
         END IF
         IF (ASSOCIATED(resp_env%rhs)) THEN
            DEALLOCATE (resp_env%rhs)
         END IF
         IF (ASSOCIATED(resp_env%sum_vpot)) THEN
            DEALLOCATE (resp_env%sum_vpot)
         END IF
         IF (ASSOCIATED(resp_env%fitpoints)) THEN
            DEALLOCATE (resp_env%fitpoints)
         END IF
         IF (ASSOCIATED(resp_env%rmin_kind)) THEN
            DEALLOCATE (resp_env%rmin_kind)
         END IF
         IF (ASSOCIATED(resp_env%rmax_kind)) THEN
            DEALLOCATE (resp_env%rmax_kind)
         END IF
         DEALLOCATE (resp_env)
      END IF
      IF (ASSOCIATED(rep_sys)) THEN
         DO i = 1, SIZE(rep_sys)
            DEALLOCATE (rep_sys(i)%p_resp%atom_surf_list)
            DEALLOCATE (rep_sys(i)%p_resp)
         END DO
         DEALLOCATE (rep_sys)
      END IF

   END SUBROUTINE resp_dealloc

! **************************************************************************************************
!> \brief initializes the resp fit. Getting the parameters
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
!> \param subsys ...
!> \param atomic_kind_set ...
!> \param cell parameters related to the simulation cell
!> \param resp_section resp section
!> \param cons_section constraints section, part of resp section
!> \param rest_section restraints section, part of resp section
! **************************************************************************************************
   SUBROUTINE init_resp(resp_env, rep_sys, subsys, atomic_kind_set, &
                        cell, resp_section, cons_section, rest_section)

      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(section_vals_type), POINTER                   :: resp_section, cons_section, rest_section

      CHARACTER(len=*), PARAMETER                        :: routineN = 'init_resp'

      INTEGER                                            :: handle, i, nrep
      INTEGER, DIMENSION(:), POINTER                     :: atom_list_cons, my_stride
      LOGICAL                                            :: explicit
      TYPE(section_vals_type), POINTER                   :: slab_section, sphere_section

      CALL timeset(routineN, handle)

      NULLIFY (atom_list_cons, my_stride, sphere_section, slab_section)

      ! get the subsections
      sphere_section => section_vals_get_subs_vals(resp_section, "SPHERE_SAMPLING")
      slab_section => section_vals_get_subs_vals(resp_section, "SLAB_SAMPLING")
      cons_section => section_vals_get_subs_vals(resp_section, "CONSTRAINT")
      rest_section => section_vals_get_subs_vals(resp_section, "RESTRAINT")

      ! get the general keywords
      CALL section_vals_val_get(resp_section, "INTEGER_TOTAL_CHARGE", &
                                l_val=resp_env%itc)
      IF (resp_env%itc) resp_env%ncons = resp_env%ncons + 1

      CALL section_vals_val_get(resp_section, "RESTRAIN_HEAVIES_TO_ZERO", &
                                l_val=resp_env%rheavies)
      IF (resp_env%rheavies) THEN
         CALL section_vals_val_get(resp_section, "RESTRAIN_HEAVIES_STRENGTH", &
                                   r_val=resp_env%rheavies_strength)
      END IF
      CALL section_vals_val_get(resp_section, "STRIDE", i_vals=my_stride)
      IF (SIZE(my_stride) /= 1 .AND. SIZE(my_stride) /= 3) &
         CALL cp_abort(__LOCATION__, "STRIDE keyword can accept only 1 (the same for X,Y,Z) "// &
                       "or 3 values. Correct your input file.")
      IF (SIZE(my_stride) == 1) THEN
         DO i = 1, 3
            resp_env%stride(i) = my_stride(1)
         END DO
      ELSE
         resp_env%stride = my_stride(1:3)
      END IF
      CALL section_vals_val_get(resp_section, "WIDTH", r_val=resp_env%eta)

      ! get if the user wants to use REPEAT method
      CALL section_vals_val_get(resp_section, "USE_REPEAT_METHOD", &
                                l_val=resp_env%use_repeat_method)
      IF (resp_env%use_repeat_method) THEN
         resp_env%ncons = resp_env%ncons + 1
         ! restrain heavies should be off
         resp_env%rheavies = .FALSE.
      END IF

      ! get and set the parameters for molecular (non-surface) systems
      ! this must come after the repeat settings being set
      CALL get_parameter_molecular_sys(resp_env, sphere_section, cell, &
                                       atomic_kind_set)

      ! get the parameter for periodic/surface systems
      CALL section_vals_get(slab_section, explicit=explicit, n_repetition=nrep)
      IF (explicit) THEN
         IF (resp_env%molecular_sys) THEN
            CALL cp_abort(__LOCATION__, &
                          "You can only use either SPHERE_SAMPLING or SLAB_SAMPLING, but "// &
                          "not both.")
         END IF
         ALLOCATE (rep_sys(nrep))
         DO i = 1, nrep
            ALLOCATE (rep_sys(i)%p_resp)
            NULLIFY (rep_sys(i)%p_resp%range_surf, rep_sys(i)%p_resp%atom_surf_list)
            CALL section_vals_val_get(slab_section, "RANGE", r_vals=rep_sys(i)%p_resp%range_surf, &
                                      i_rep_section=i)
            CALL section_vals_val_get(slab_section, "LENGTH", r_val=rep_sys(i)%p_resp%length, &
                                      i_rep_section=i)
            CALL section_vals_val_get(slab_section, "SURF_DIRECTION", &
                                      i_rep_section=i, i_val=rep_sys(i)%p_resp%my_fit)
            IF (ANY(rep_sys(i)%p_resp%range_surf < 0.0_dp)) THEN
               CPABORT("Numbers in RANGE in SLAB_SAMPLING cannot be negative.")
            END IF
            IF (rep_sys(i)%p_resp%length <= EPSILON(0.0_dp)) THEN
               CPABORT("Parameter LENGTH in SLAB_SAMPLING has to be larger than zero.")
            END IF
            !list of atoms specifying the surface
            CALL build_atom_list(slab_section, subsys, rep_sys(i)%p_resp%atom_surf_list, rep=i)
         END DO
      END IF

      ! get the parameters for the constraint and restraint sections
      CALL section_vals_get(cons_section, explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_get(cons_section, n_repetition=resp_env%ncons_sec)
         DO i = 1, resp_env%ncons_sec
            CALL section_vals_val_get(cons_section, "EQUAL_CHARGES", &
                                      l_val=resp_env%equal_charges, explicit=explicit)
            IF (.NOT. explicit) CYCLE
            CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
            !instead of using EQUAL_CHARGES the constraint sections could be repeated
            resp_env%ncons = resp_env%ncons + SIZE(atom_list_cons) - 2
            DEALLOCATE (atom_list_cons)
         END DO
      END IF
      CALL section_vals_get(rest_section, explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_get(rest_section, n_repetition=resp_env%nrest_sec)
      END IF
      resp_env%ncons = resp_env%ncons + resp_env%ncons_sec
      resp_env%nres = resp_env%nres + resp_env%nrest_sec

      CALL timestop(handle)

   END SUBROUTINE init_resp

! **************************************************************************************************
!> \brief getting the parameters for nonperiodic/non-surface systems
!> \param resp_env the resp environment
!> \param sphere_section input section setting parameters for sampling
!>        fitting in spheres around the atom
!> \param cell parameters related to the simulation cell
!> \param atomic_kind_set ...
! **************************************************************************************************
   SUBROUTINE get_parameter_molecular_sys(resp_env, sphere_section, cell, &
                                          atomic_kind_set)

      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(section_vals_type), POINTER                   :: sphere_section
      TYPE(cell_type), POINTER                           :: cell
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set

      CHARACTER(LEN=2)                                   :: symbol
      CHARACTER(LEN=default_string_length)               :: missing_rmax, missing_rmin
      CHARACTER(LEN=default_string_length), &
         DIMENSION(:), POINTER                           :: tmpstringlist
      INTEGER                                            :: ikind, j, kind_number, n_rmax_missing, &
                                                            n_rmin_missing, nkind, nrep_rmax, &
                                                            nrep_rmin, z
      LOGICAL                                            :: explicit, has_rmax, has_rmin
      LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: rmax_is_set, rmin_is_set
      REAL(KIND=dp)                                      :: auto_rmax_scale, auto_rmin_scale, rmax, &
                                                            rmin
      REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
      TYPE(atomic_kind_type), POINTER                    :: atomic_kind

      nrep_rmin = 0
      nrep_rmax = 0
      nkind = SIZE(atomic_kind_set)

      has_rmin = .FALSE.
      has_rmax = .FALSE.

      CALL section_vals_get(sphere_section, explicit=explicit)
      IF (explicit) THEN
         resp_env%molecular_sys = .TRUE.
         CALL section_vals_val_get(sphere_section, "AUTO_VDW_RADII_TABLE", &
                                   i_val=resp_env%auto_vdw_radii_table)
         CALL section_vals_val_get(sphere_section, "AUTO_RMIN_SCALE", r_val=auto_rmin_scale)
         CALL section_vals_val_get(sphere_section, "AUTO_RMAX_SCALE", r_val=auto_rmax_scale)
         CALL section_vals_val_get(sphere_section, "RMIN", explicit=has_rmin, r_val=rmin)
         CALL section_vals_val_get(sphere_section, "RMAX", explicit=has_rmax, r_val=rmax)
         CALL section_vals_val_get(sphere_section, "RMIN_KIND", n_rep_val=nrep_rmin)
         CALL section_vals_val_get(sphere_section, "RMAX_KIND", n_rep_val=nrep_rmax)
         ALLOCATE (resp_env%rmin_kind(nkind))
         ALLOCATE (resp_env%rmax_kind(nkind))
         resp_env%rmin_kind = 0.0_dp
         resp_env%rmax_kind = 0.0_dp
         ALLOCATE (rmin_is_set(nkind))
         ALLOCATE (rmax_is_set(nkind))
         rmin_is_set = .FALSE.
         rmax_is_set = .FALSE.
         ! define rmin_kind and rmax_kind to predefined vdW radii
         DO ikind = 1, nkind
            atomic_kind => atomic_kind_set(ikind)
            CALL get_atomic_kind(atomic_kind, &
                                 element_symbol=symbol, &
                                 kind_number=kind_number, &
                                 z=z)
            SELECT CASE (resp_env%auto_vdw_radii_table)
            CASE (use_cambridge_vdw_radii)
               CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
               rmin_is_set(kind_number) = .TRUE.
            CASE (use_uff_vdw_radii)
               CALL cite_reference(Rappe1992)
               CALL get_uff_vdw_radius(z, radius=resp_env%rmin_kind(kind_number), &
                                       found=rmin_is_set(kind_number))
            CASE DEFAULT
               CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
               rmin_is_set(kind_number) = .TRUE.
            END SELECT
            IF (rmin_is_set(kind_number)) THEN
               resp_env%rmin_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
                                                                 "angstrom")
               resp_env%rmin_kind(kind_number) = resp_env%rmin_kind(kind_number)*auto_rmin_scale
               ! set RMAX_KIND accourding by scaling RMIN_KIND
               resp_env%rmax_kind(kind_number) = &
                  MAX(resp_env%rmin_kind(kind_number), &
                      resp_env%rmin_kind(kind_number)*auto_rmax_scale)
               rmax_is_set(kind_number) = .TRUE.
            END IF
         END DO
         ! if RMIN or RMAX are present, overwrite the rmin_kind(:) and
         ! rmax_kind(:) to those values
         IF (has_rmin) THEN
            resp_env%rmin_kind = rmin
            rmin_is_set = .TRUE.
         END IF
         IF (has_rmax) THEN
            resp_env%rmax_kind = rmax
            rmax_is_set = .TRUE.
         END IF
         ! if RMIN_KIND's or RMAX_KIND's are present, overwrite the
         ! rmin_kinds(:) or rmax_kind(:) to those values
         DO j = 1, nrep_rmin
            CALL section_vals_val_get(sphere_section, "RMIN_KIND", i_rep_val=j, &
                                      c_vals=tmpstringlist)
            DO ikind = 1, nkind
               atomic_kind => atomic_kind_set(ikind)
               CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
               IF (TRIM(tmpstringlist(2)) == TRIM(symbol)) THEN
                  READ (tmpstringlist(1), *) resp_env%rmin_kind(kind_number)
                  resp_env%rmin_kind(kind_number) = &
                     cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
                                     "angstrom")
                  rmin_is_set(kind_number) = .TRUE.
               END IF
            END DO
         END DO
         DO j = 1, nrep_rmax
            CALL section_vals_val_get(sphere_section, "RMAX_KIND", i_rep_val=j, &
                                      c_vals=tmpstringlist)
            DO ikind = 1, nkind
               atomic_kind => atomic_kind_set(ikind)
               CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
               IF (TRIM(tmpstringlist(2)) == TRIM(symbol)) THEN
                  READ (tmpstringlist(1), *) resp_env%rmax_kind(kind_number)
                  resp_env%rmax_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmax_kind(kind_number), &
                                                                    "angstrom")
                  rmax_is_set(kind_number) = .TRUE.
               END IF
            END DO
         END DO
         ! check if rmin and rmax are set for each kind
         n_rmin_missing = 0
         n_rmax_missing = 0
         missing_rmin = ""
         missing_rmax = ""
         DO ikind = 1, nkind
            atomic_kind => atomic_kind_set(ikind)
            CALL get_atomic_kind(atomic_kind, &
                                 element_symbol=symbol, &
                                 kind_number=kind_number)
            IF (.NOT. rmin_is_set(kind_number)) THEN
               n_rmin_missing = n_rmin_missing + 1
               missing_rmin = TRIM(missing_rmin)//" "//TRIM(symbol)//","
            END IF
            IF (.NOT. rmax_is_set(kind_number)) THEN
               n_rmax_missing = n_rmax_missing + 1
               missing_rmax = TRIM(missing_rmax)//" "//TRIM(symbol)//","
            END IF
         END DO
         IF (n_rmin_missing > 0) THEN
            CALL cp_warn(__LOCATION__, &
                         "RMIN for the following elements are missing: "// &
                         TRIM(missing_rmin)// &
                         " please set these values manually using "// &
                         "RMIN_KIND in SPHERE_SAMPLING section")
         END IF
         IF (n_rmax_missing > 0) THEN
            CALL cp_warn(__LOCATION__, &
                         "RMAX for the following elements are missing: "// &
                         TRIM(missing_rmax)// &
                         " please set these values manually using "// &
                         "RMAX_KIND in SPHERE_SAMPLING section")
         END IF
         IF (n_rmin_missing > 0 .OR. &
             n_rmax_missing > 0) THEN
            CPABORT("Insufficient data for RMIN or RMAX")
         END IF

         CALL get_cell(cell=cell, h=hmat)
         resp_env%box_hi = (/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
         resp_env%box_low = 0.0_dp
         CALL section_vals_val_get(sphere_section, "X_HI", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "X_HI", &
                                                 r_val=resp_env%box_hi(1))
         CALL section_vals_val_get(sphere_section, "X_LOW", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "X_LOW", &
                                                 r_val=resp_env%box_low(1))
         CALL section_vals_val_get(sphere_section, "Y_HI", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "Y_HI", &
                                                 r_val=resp_env%box_hi(2))
         CALL section_vals_val_get(sphere_section, "Y_LOW", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "Y_LOW", &
                                                 r_val=resp_env%box_low(2))
         CALL section_vals_val_get(sphere_section, "Z_HI", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "Z_HI", &
                                                 r_val=resp_env%box_hi(3))
         CALL section_vals_val_get(sphere_section, "Z_LOW", explicit=explicit)
         IF (explicit) CALL section_vals_val_get(sphere_section, "Z_LOW", &
                                                 r_val=resp_env%box_low(3))

         DEALLOCATE (rmin_is_set)
         DEALLOCATE (rmax_is_set)
      END IF

   END SUBROUTINE get_parameter_molecular_sys

! **************************************************************************************************
!> \brief building atom lists for different sections of RESP
!> \param section input section
!> \param subsys ...
!> \param atom_list list of atoms for restraints, constraints and fit point
!>        sampling for slab-like systems
!> \param rep input section can be repeated, this param defines for which
!>        repetition of the input section the atom_list is built
! **************************************************************************************************
   SUBROUTINE build_atom_list(section, subsys, atom_list, rep)

      TYPE(section_vals_type), POINTER                   :: section
      TYPE(qs_subsys_type), POINTER                      :: subsys
      INTEGER, DIMENSION(:), POINTER                     :: atom_list
      INTEGER, INTENT(IN), OPTIONAL                      :: rep

      CHARACTER(len=*), PARAMETER                        :: routineN = 'build_atom_list'

      INTEGER                                            :: atom_a, atom_b, handle, i, irep, j, &
                                                            max_index, n_var, num_atom
      INTEGER, DIMENSION(:), POINTER                     :: indexes
      LOGICAL                                            :: index_in_range

      CALL timeset(routineN, handle)

      NULLIFY (indexes)
      irep = 1
      IF (PRESENT(rep)) irep = rep

      CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
                                n_rep_val=n_var)
      num_atom = 0
      DO i = 1, n_var
         CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
                                   i_rep_val=i, i_vals=indexes)
         num_atom = num_atom + SIZE(indexes)
      END DO
      ALLOCATE (atom_list(num_atom))
      atom_list = 0
      num_atom = 1
      DO i = 1, n_var
         CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
                                   i_rep_val=i, i_vals=indexes)
         atom_list(num_atom:num_atom + SIZE(indexes) - 1) = indexes(:)
         num_atom = num_atom + SIZE(indexes)
      END DO
      !check atom list
      num_atom = num_atom - 1
      CALL qs_subsys_get(subsys, nparticle=max_index)
      CPASSERT(SIZE(atom_list) /= 0)
      index_in_range = (MAXVAL(atom_list) <= max_index) &
                       .AND. (MINVAL(atom_list) > 0)
      CPASSERT(index_in_range)
      DO i = 1, num_atom
         DO j = i + 1, num_atom
            atom_a = atom_list(i)
            atom_b = atom_list(j)
            IF (atom_a == atom_b) &
               CPABORT("There are atoms doubled in atom list for RESP.")
         END DO
      END DO

      CALL timestop(handle)

   END SUBROUTINE build_atom_list

! **************************************************************************************************
!> \brief build matrix and vector for nonperiodic RESP fitting
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param atomic_kind_set ...
!> \param particles ...
!> \param cell parameters related to the simulation cell
!> \param matrix coefficient matrix of the linear system of equations
!> \param rhs vector of the linear system of equations
!> \param natom number of atoms
! **************************************************************************************************
   SUBROUTINE calc_resp_matrix_nonper(qs_env, resp_env, atomic_kind_set, particles, &
                                      cell, matrix, rhs, natom)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(cell_type), POINTER                           :: cell
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: matrix
      REAL(KIND=dp), DIMENSION(:), POINTER               :: rhs
      INTEGER, INTENT(IN)                                :: natom

      CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_matrix_nonper'

      INTEGER                                            :: bo(2, 3), gbo(2, 3), handle, i, ikind, &
                                                            jx, jy, jz, k, kind_number, l, m, &
                                                            nkind, now, np(3), p
      LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: not_in_range
      REAL(KIND=dp)                                      :: delta, dh(3, 3), dvol, r(3), rmax, rmin, &
                                                            vec(3), vec_pbc(3), vj
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dist
      REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat, hmat_inv
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_pw

      CALL timeset(routineN, handle)

      NULLIFY (particle_set, v_hartree_pw)
      delta = 1.0E-13_dp

      CALL get_cell(cell=cell, h=hmat, h_inv=hmat_inv)

      IF (.NOT. cell%orthorhombic) THEN
         CALL cp_abort(__LOCATION__, &
                       "Nonperiodic solution for RESP charges only"// &
                       " implemented for orthorhombic cells!")
      END IF
      IF (.NOT. resp_env%molecular_sys) THEN
         CALL cp_abort(__LOCATION__, &
                       "Nonperiodic solution for RESP charges (i.e. nonperiodic"// &
                       " Poisson solver) can only be used with section SPHERE_SAMPLING")
      END IF
      IF (resp_env%use_repeat_method) THEN
         CALL cp_abort(__LOCATION__, &
                       "REPEAT method only reasonable for periodic RESP fitting")
      END IF
      CALL get_qs_env(qs_env, particle_set=particle_set, v_hartree_rspace=v_hartree_pw)

      bo = v_hartree_pw%pw_grid%bounds_local
      gbo = v_hartree_pw%pw_grid%bounds
      np = v_hartree_pw%pw_grid%npts
      dh = v_hartree_pw%pw_grid%dh
      dvol = v_hartree_pw%pw_grid%dvol
      nkind = SIZE(atomic_kind_set)

      ALLOCATE (dist(natom))
      ALLOCATE (not_in_range(natom, 2))

      ! store fitting points to calculate the RMS and RRMS later
      IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
         now = 1000
         ALLOCATE (resp_env%fitpoints(3, now))
      ELSE
         now = SIZE(resp_env%fitpoints, 2)
      END IF

      DO jz = bo(1, 3), bo(2, 3)
      DO jy = bo(1, 2), bo(2, 2)
      DO jx = bo(1, 1), bo(2, 1)
         IF (.NOT. (MODULO(jz, resp_env%stride(3)) == 0)) CYCLE
         IF (.NOT. (MODULO(jy, resp_env%stride(2)) == 0)) CYCLE
         IF (.NOT. (MODULO(jx, resp_env%stride(1)) == 0)) CYCLE
         !bounds bo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
         l = jx - gbo(1, 1)
         k = jy - gbo(1, 2)
         p = jz - gbo(1, 3)
         r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
         r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
         r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
         IF (r(3) < resp_env%box_low(3) .OR. r(3) > resp_env%box_hi(3)) CYCLE
         IF (r(2) < resp_env%box_low(2) .OR. r(2) > resp_env%box_hi(2)) CYCLE
         IF (r(1) < resp_env%box_low(1) .OR. r(1) > resp_env%box_hi(1)) CYCLE
         ! compute distance from the grid point to all atoms
         not_in_range = .FALSE.
         DO i = 1, natom
            vec = r - particles%els(i)%r
            vec_pbc(1) = vec(1) - hmat(1, 1)*ANINT(hmat_inv(1, 1)*vec(1))
            vec_pbc(2) = vec(2) - hmat(2, 2)*ANINT(hmat_inv(2, 2)*vec(2))
            vec_pbc(3) = vec(3) - hmat(3, 3)*ANINT(hmat_inv(3, 3)*vec(3))
            dist(i) = SQRT(SUM(vec_pbc**2))
            CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
                                 kind_number=kind_number)
            DO ikind = 1, nkind
               IF (ikind == kind_number) THEN
                  rmin = resp_env%rmin_kind(ikind)
                  rmax = resp_env%rmax_kind(ikind)
                  EXIT
               END IF
            END DO
            IF (dist(i) < rmin + delta) not_in_range(i, 1) = .TRUE.
            IF (dist(i) > rmax - delta) not_in_range(i, 2) = .TRUE.
         END DO
         ! check if the point is sufficiently close and  far. if OK, we can use
         ! the point for fitting, add/subtract 1.0E-13 to get rid of rounding errors when shifting atoms
         IF (ANY(not_in_range(:, 1)) .OR. ALL(not_in_range(:, 2))) CYCLE
         resp_env%npoints_proc = resp_env%npoints_proc + 1
         IF (resp_env%npoints_proc > now) THEN
            now = 2*now
            CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
         END IF
         resp_env%fitpoints(1, resp_env%npoints_proc) = jx
         resp_env%fitpoints(2, resp_env%npoints_proc) = jy
         resp_env%fitpoints(3, resp_env%npoints_proc) = jz
         ! correct for the fact that v_hartree is scaled by dvol, and has the opposite sign
         IF (qs_env%qmmm) THEN
            ! If it's a QM/MM run let's remove the contribution of the MM potential out of the Hartree pot
            vj = -v_hartree_pw%array(jx, jy, jz)/dvol + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
         ELSE
            vj = -v_hartree_pw%array(jx, jy, jz)/dvol
         END IF
         dist(:) = 1.0_dp/dist(:)

         DO i = 1, natom
            DO m = 1, natom
               matrix(m, i) = matrix(m, i) + 2.0_dp*dist(i)*dist(m)
            END DO
            rhs(i) = rhs(i) + 2.0_dp*vj*dist(i)
         END DO
      END DO
      END DO
      END DO

      resp_env%npoints = resp_env%npoints_proc
      CALL v_hartree_pw%pw_grid%para%group%sum(resp_env%npoints)
      CALL v_hartree_pw%pw_grid%para%group%sum(matrix)
      CALL v_hartree_pw%pw_grid%para%group%sum(rhs)
      !weighted sum
      matrix = matrix/resp_env%npoints
      rhs = rhs/resp_env%npoints

      DEALLOCATE (dist)
      DEALLOCATE (not_in_range)

      CALL timestop(handle)

   END SUBROUTINE calc_resp_matrix_nonper

! **************************************************************************************************
!> \brief build matrix and vector for periodic RESP fitting
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
!> \param particles ...
!> \param cell parameters related to the simulation cell
!> \param natom number of atoms
! **************************************************************************************************
   SUBROUTINE calc_resp_matrix_periodic(qs_env, resp_env, rep_sys, particles, cell, &
                                        natom)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(cell_type), POINTER                           :: cell
      INTEGER, INTENT(IN)                                :: natom

      CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_matrix_periodic'

      INTEGER                                            :: handle, i, ip, j, jx, jy, jz
      INTEGER, DIMENSION(3)                              :: periodic
      REAL(KIND=dp)                                      :: normalize_factor
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: vpot
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(pw_c1d_gs_type)                               :: rho_ga, va_gspace
      TYPE(pw_env_type), POINTER                         :: pw_env
      TYPE(pw_poisson_type), POINTER                     :: poisson_env
      TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      TYPE(pw_r3d_rs_type)                               :: va_rspace

      CALL timeset(routineN, handle)

      NULLIFY (pw_env, para_env, auxbas_pw_pool, poisson_env)

      CALL get_cell(cell=cell, periodic=periodic)

      IF (.NOT. ALL(periodic /= 0)) THEN
         CALL cp_abort(__LOCATION__, &
                       "Periodic solution for RESP (with periodic Poisson solver)"// &
                       " can only be obtained with a cell that has XYZ periodicity")
      END IF

      CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)

      CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
                      poisson_env=poisson_env)
      CALL auxbas_pw_pool%create_pw(rho_ga)
      CALL auxbas_pw_pool%create_pw(va_gspace)
      CALL auxbas_pw_pool%create_pw(va_rspace)

      !get fitting points and store them in resp_env%fitpoints
      CALL get_fitting_points(qs_env, resp_env, rep_sys, particles=particles, &
                              cell=cell)
      ALLOCATE (vpot(resp_env%npoints_proc, natom))
      normalize_factor = SQRT((resp_env%eta/pi)**3)

      DO i = 1, natom
         !collocate gaussian for each atom
         CALL pw_zero(rho_ga)
         CALL calculate_rho_resp_single(rho_ga, qs_env, resp_env%eta, i)
         !calculate potential va and store the part needed for fitting in vpot
         CALL pw_zero(va_gspace)
         CALL pw_poisson_solve(poisson_env, rho_ga, vhartree=va_gspace)
         CALL pw_zero(va_rspace)
         CALL pw_transfer(va_gspace, va_rspace)
         CALL pw_scale(va_rspace, normalize_factor)
         DO ip = 1, resp_env%npoints_proc
            jx = resp_env%fitpoints(1, ip)
            jy = resp_env%fitpoints(2, ip)
            jz = resp_env%fitpoints(3, ip)
            vpot(ip, i) = va_rspace%array(jx, jy, jz)
         END DO
      END DO

      CALL va_gspace%release()
      CALL va_rspace%release()
      CALL rho_ga%release()

      DO i = 1, natom
         DO j = 1, natom
            ! calculate matrix
            resp_env%matrix(i, j) = resp_env%matrix(i, j) + 2.0_dp*SUM(vpot(:, i)*vpot(:, j))
         END DO
         ! calculate vector resp_env%rhs
         CALL calculate_rhs(qs_env, resp_env, resp_env%rhs(i), vpot(:, i))
      END DO

      CALL para_env%sum(resp_env%matrix)
      CALL para_env%sum(resp_env%rhs)
      !weighted sum
      resp_env%matrix = resp_env%matrix/resp_env%npoints
      resp_env%rhs = resp_env%rhs/resp_env%npoints

      ! REPEAT stuff
      IF (resp_env%use_repeat_method) THEN
         ! sum over selected points of single Gaussian potential vpot
         DO i = 1, natom
            resp_env%sum_vpot(i) = 2.0_dp*accurate_sum(vpot(:, i))/resp_env%npoints
         END DO
         CALL para_env%sum(resp_env%sum_vpot)
         CALL para_env%sum(resp_env%sum_vhartree)
         resp_env%sum_vhartree = 2.0_dp*resp_env%sum_vhartree/resp_env%npoints
      END IF

      DEALLOCATE (vpot)
      CALL timestop(handle)

   END SUBROUTINE calc_resp_matrix_periodic

! **************************************************************************************************
!> \brief get RESP fitting points for the periodic fitting
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
!> \param particles ...
!> \param cell parameters related to the simulation cell
! **************************************************************************************************
   SUBROUTINE get_fitting_points(qs_env, resp_env, rep_sys, particles, cell)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(cell_type), POINTER                           :: cell

      CHARACTER(len=*), PARAMETER :: routineN = 'get_fitting_points'

      INTEGER                                            :: bo(2, 3), gbo(2, 3), handle, i, iatom, &
                                                            ikind, in_x, in_y, in_z, jx, jy, jz, &
                                                            k, kind_number, l, m, natom, nkind, &
                                                            now, p
      LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: not_in_range
      REAL(KIND=dp)                                      :: delta, dh(3, 3), r(3), rmax, rmin, &
                                                            vec_pbc(3)
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dist
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_pw

      CALL timeset(routineN, handle)

      NULLIFY (atomic_kind_set, v_hartree_pw, para_env, particle_set)
      delta = 1.0E-13_dp

      CALL get_qs_env(qs_env, &
                      particle_set=particle_set, &
                      atomic_kind_set=atomic_kind_set, &
                      para_env=para_env, &
                      v_hartree_rspace=v_hartree_pw)

      bo = v_hartree_pw%pw_grid%bounds_local
      gbo = v_hartree_pw%pw_grid%bounds
      dh = v_hartree_pw%pw_grid%dh
      natom = SIZE(particles%els)
      nkind = SIZE(atomic_kind_set)

      IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
         now = 1000
         ALLOCATE (resp_env%fitpoints(3, now))
      ELSE
         now = SIZE(resp_env%fitpoints, 2)
      END IF

      ALLOCATE (dist(natom))
      ALLOCATE (not_in_range(natom, 2))

      !every proc gets another bo, grid is distributed
      DO jz = bo(1, 3), bo(2, 3)
         IF (.NOT. (MODULO(jz, resp_env%stride(3)) == 0)) CYCLE
         DO jy = bo(1, 2), bo(2, 2)
            IF (.NOT. (MODULO(jy, resp_env%stride(2)) == 0)) CYCLE
            DO jx = bo(1, 1), bo(2, 1)
               IF (.NOT. (MODULO(jx, resp_env%stride(1)) == 0)) CYCLE
               !bounds gbo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
               l = jx - gbo(1, 1)
               k = jy - gbo(1, 2)
               p = jz - gbo(1, 3)
               r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
               r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
               r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
               IF (resp_env%molecular_sys) THEN
                  not_in_range = .FALSE.
                  DO m = 1, natom
                     vec_pbc = pbc(r, particles%els(m)%r, cell)
                     dist(m) = SQRT(SUM(vec_pbc**2))
                     CALL get_atomic_kind(atomic_kind=particle_set(m)%atomic_kind, &
                                          kind_number=kind_number)
                     DO ikind = 1, nkind
                        IF (ikind == kind_number) THEN
                           rmin = resp_env%rmin_kind(ikind)
                           rmax = resp_env%rmax_kind(ikind)
                           EXIT
                        END IF
                     END DO
                     IF (dist(m) < rmin + delta) not_in_range(m, 1) = .TRUE.
                     IF (dist(m) > rmax - delta) not_in_range(m, 2) = .TRUE.
                  END DO
                  IF (ANY(not_in_range(:, 1)) .OR. ALL(not_in_range(:, 2))) CYCLE
               ELSE
                  DO i = 1, SIZE(rep_sys)
                     DO m = 1, SIZE(rep_sys(i)%p_resp%atom_surf_list)
                        in_z = 0
                        in_y = 0
                        in_x = 0
                        iatom = rep_sys(i)%p_resp%atom_surf_list(m)
                        SELECT CASE (rep_sys(i)%p_resp%my_fit)
                        CASE (do_resp_x_dir, do_resp_y_dir, do_resp_z_dir)
                           vec_pbc = pbc(particles%els(iatom)%r, r, cell)
                        CASE (do_resp_minus_x_dir, do_resp_minus_y_dir, do_resp_minus_z_dir)
                           vec_pbc = pbc(r, particles%els(iatom)%r, cell)
                        END SELECT
                        SELECT CASE (rep_sys(i)%p_resp%my_fit)
                           !subtract delta=1.0E-13 to get rid of rounding errors when shifting atoms
                        CASE (do_resp_x_dir, do_resp_minus_x_dir)
                           IF (ABS(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
                           IF (ABS(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
                           IF (vec_pbc(1) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
                               vec_pbc(1) < rep_sys(i)%p_resp%range_surf(2) - delta) in_x = 1
                        CASE (do_resp_y_dir, do_resp_minus_y_dir)
                           IF (ABS(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
                           IF (vec_pbc(2) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
                               vec_pbc(2) < rep_sys(i)%p_resp%range_surf(2) - delta) in_y = 1
                           IF (ABS(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
                        CASE (do_resp_z_dir, do_resp_minus_z_dir)
                           IF (vec_pbc(3) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
                               vec_pbc(3) < rep_sys(i)%p_resp%range_surf(2) - delta) in_z = 1
                           IF (ABS(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
                           IF (ABS(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
                        END SELECT
                        IF (in_z*in_y*in_x == 1) EXIT
                     END DO
                     IF (in_z*in_y*in_x == 1) EXIT
                  END DO
                  IF (in_z*in_y*in_x == 0) CYCLE
               END IF
               resp_env%npoints_proc = resp_env%npoints_proc + 1
               IF (resp_env%npoints_proc > now) THEN
                  now = 2*now
                  CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
               END IF
               resp_env%fitpoints(1, resp_env%npoints_proc) = jx
               resp_env%fitpoints(2, resp_env%npoints_proc) = jy
               resp_env%fitpoints(3, resp_env%npoints_proc) = jz
            END DO
         END DO
      END DO

      resp_env%npoints = resp_env%npoints_proc
      CALL para_env%sum(resp_env%npoints)

      DEALLOCATE (dist)
      DEALLOCATE (not_in_range)

      CALL timestop(handle)

   END SUBROUTINE get_fitting_points

! **************************************************************************************************
!> \brief calculate vector rhs
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param rhs vector
!> \param vpot single gaussian potential
! **************************************************************************************************
   SUBROUTINE calculate_rhs(qs_env, resp_env, rhs, vpot)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      REAL(KIND=dp), INTENT(INOUT)                       :: rhs
      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: vpot

      CHARACTER(len=*), PARAMETER                        :: routineN = 'calculate_rhs'

      INTEGER                                            :: handle, ip, jx, jy, jz
      REAL(KIND=dp)                                      :: dvol
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: vhartree
      TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_pw

      CALL timeset(routineN, handle)

      NULLIFY (v_hartree_pw)
      CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_pw)
      dvol = v_hartree_pw%pw_grid%dvol
      ALLOCATE (vhartree(resp_env%npoints_proc))
      vhartree = 0.0_dp

      !multiply v_hartree and va_rspace and calculate the vector rhs
      !taking into account that v_hartree has opposite site; remove v_qmmm
      DO ip = 1, resp_env%npoints_proc
         jx = resp_env%fitpoints(1, ip)
         jy = resp_env%fitpoints(2, ip)
         jz = resp_env%fitpoints(3, ip)
         vhartree(ip) = -v_hartree_pw%array(jx, jy, jz)/dvol
         IF (qs_env%qmmm) THEN
            !taking into account that v_qmmm has also opposite sign
            vhartree(ip) = vhartree(ip) + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
         END IF
         rhs = rhs + 2.0_dp*vhartree(ip)*vpot(ip)
      END DO

      IF (resp_env%use_repeat_method) THEN
         resp_env%sum_vhartree = accurate_sum(vhartree)
      END IF

      DEALLOCATE (vhartree)

      CALL timestop(handle)

   END SUBROUTINE calculate_rhs

! **************************************************************************************************
!> \brief print the atom coordinates and the coordinates of the fitting points
!>        to an xyz file
!> \param qs_env the qs environment
!> \param resp_env the resp environment
! **************************************************************************************************
   SUBROUTINE print_fitting_points(qs_env, resp_env)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env

      CHARACTER(len=*), PARAMETER :: routineN = 'print_fitting_points'

      CHARACTER(LEN=2)                                   :: element_symbol
      CHARACTER(LEN=default_path_length)                 :: filename
      INTEGER                                            :: gbo(2, 3), handle, i, iatom, ip, jx, jy, &
                                                            jz, k, l, my_pos, nobjects, &
                                                            output_unit, p
      INTEGER, DIMENSION(:), POINTER                     :: tmp_npoints, tmp_size
      INTEGER, DIMENSION(:, :), POINTER                  :: tmp_points
      REAL(KIND=dp)                                      :: conv, dh(3, 3), r(3)
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(mp_request_type), DIMENSION(6)                :: req
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_pw
      TYPE(section_vals_type), POINTER                   :: input, print_key, resp_section

      CALL timeset(routineN, handle)

      NULLIFY (para_env, input, logger, resp_section, print_key, particle_set, tmp_size, &
               tmp_points, tmp_npoints, v_hartree_pw)

      CALL get_qs_env(qs_env, input=input, para_env=para_env, &
                      particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
      conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
      gbo = v_hartree_pw%pw_grid%bounds
      dh = v_hartree_pw%pw_grid%dh
      nobjects = SIZE(particle_set) + resp_env%npoints

      resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
      print_key => section_vals_get_subs_vals(resp_section, "PRINT%COORD_FIT_POINTS")
      logger => cp_get_default_logger()
      output_unit = cp_print_key_unit_nr(logger, resp_section, &
                                         "PRINT%COORD_FIT_POINTS", &
                                         extension=".xyz", &
                                         file_status="REPLACE", &
                                         file_action="WRITE", &
                                         file_form="FORMATTED")

      IF (BTEST(cp_print_key_should_output(logger%iter_info, &
                                           resp_section, "PRINT%COORD_FIT_POINTS"), &
                cp_p_file)) THEN
         IF (output_unit > 0) THEN
            filename = cp_print_key_generate_filename(logger, &
                                                      print_key, extension=".xyz", &
                                                      my_local=.FALSE.)
            WRITE (unit=output_unit, FMT="(I12,/)") nobjects
            DO iatom = 1, SIZE(particle_set)
               CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
                                    element_symbol=element_symbol)
               WRITE (UNIT=output_unit, FMT="(A,1X,3F10.5)") element_symbol, &
                  particle_set(iatom)%r(1:3)*conv
            END DO
            !printing points of proc which is doing the output (should be proc 0)
            DO ip = 1, resp_env%npoints_proc
               jx = resp_env%fitpoints(1, ip)
               jy = resp_env%fitpoints(2, ip)
               jz = resp_env%fitpoints(3, ip)
               l = jx - gbo(1, 1)
               k = jy - gbo(1, 2)
               p = jz - gbo(1, 3)
               r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
               r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
               r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
               r(:) = r(:)*conv
               WRITE (UNIT=output_unit, FMT="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
            END DO
         END IF

         ALLOCATE (tmp_size(1))
         ALLOCATE (tmp_npoints(1))

         !sending data of all other procs to proc which makes the output (proc 0)
         IF (output_unit > 0) THEN
            my_pos = para_env%mepos
            DO i = 1, para_env%num_pe
               IF (my_pos == i - 1) CYCLE
               CALL para_env%irecv(msgout=tmp_size, source=i - 1, &
                                   request=req(1))
               CALL req(1)%wait()
               ALLOCATE (tmp_points(3, tmp_size(1)))
               CALL para_env%irecv(msgout=tmp_points, source=i - 1, &
                                   request=req(3))
               CALL req(3)%wait()
               CALL para_env%irecv(msgout=tmp_npoints, source=i - 1, &
                                   request=req(5))
               CALL req(5)%wait()
               DO ip = 1, tmp_npoints(1)
                  jx = tmp_points(1, ip)
                  jy = tmp_points(2, ip)
                  jz = tmp_points(3, ip)
                  l = jx - gbo(1, 1)
                  k = jy - gbo(1, 2)
                  p = jz - gbo(1, 3)
                  r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
                  r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
                  r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
                  r(:) = r(:)*conv
                  WRITE (UNIT=output_unit, FMT="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
               END DO
               DEALLOCATE (tmp_points)
            END DO
         ELSE
            tmp_size(1) = SIZE(resp_env%fitpoints, 2)
            !para_env%source should be 0
            CALL para_env%isend(msgin=tmp_size, dest=para_env%source, &
                                request=req(2))
            CALL req(2)%wait()
            CALL para_env%isend(msgin=resp_env%fitpoints, dest=para_env%source, &
                                request=req(4))
            CALL req(4)%wait()
            tmp_npoints(1) = resp_env%npoints_proc
            CALL para_env%isend(msgin=tmp_npoints, dest=para_env%source, &
                                request=req(6))
            CALL req(6)%wait()
         END IF

         DEALLOCATE (tmp_size)
         DEALLOCATE (tmp_npoints)
      END IF

      CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
                                        "PRINT%COORD_FIT_POINTS")

      CALL timestop(handle)

   END SUBROUTINE print_fitting_points

! **************************************************************************************************
!> \brief add restraints and constraints
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param rest_section input section for restraints
!> \param subsys ...
!> \param natom number of atoms
!> \param cons_section input section for constraints
!> \param particle_set ...
! **************************************************************************************************
   SUBROUTINE add_restraints_and_constraints(qs_env, resp_env, rest_section, &
                                             subsys, natom, cons_section, particle_set)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(section_vals_type), POINTER                   :: rest_section
      TYPE(qs_subsys_type), POINTER                      :: subsys
      INTEGER, INTENT(IN)                                :: natom
      TYPE(section_vals_type), POINTER                   :: cons_section
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set

      CHARACTER(len=*), PARAMETER :: routineN = 'add_restraints_and_constraints'

      INTEGER                                            :: handle, i, k, m, ncons_v, z
      INTEGER, DIMENSION(:), POINTER                     :: atom_list_cons, atom_list_res
      LOGICAL                                            :: explicit_coeff
      REAL(KIND=dp)                                      :: my_atom_coef(2), strength, TARGET
      REAL(KIND=dp), DIMENSION(:), POINTER               :: atom_coef
      TYPE(dft_control_type), POINTER                    :: dft_control

      CALL timeset(routineN, handle)

      NULLIFY (atom_coef, atom_list_res, atom_list_cons, dft_control)

      CALL get_qs_env(qs_env, dft_control=dft_control)

      !*** add the restraints
      DO i = 1, resp_env%nrest_sec
         CALL section_vals_val_get(rest_section, "TARGET", i_rep_section=i, r_val=TARGET)
         CALL section_vals_val_get(rest_section, "STRENGTH", i_rep_section=i, r_val=strength)
         CALL build_atom_list(rest_section, subsys, atom_list_res, i)
         CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, explicit=explicit_coeff)
         IF (explicit_coeff) THEN
            CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
            CPASSERT(SIZE(atom_list_res) == SIZE(atom_coef))
         END IF
         DO m = 1, SIZE(atom_list_res)
            IF (explicit_coeff) THEN
               DO k = 1, SIZE(atom_list_res)
                  resp_env%matrix(atom_list_res(m), atom_list_res(k)) = &
                     resp_env%matrix(atom_list_res(m), atom_list_res(k)) + &
                     atom_coef(m)*atom_coef(k)*2.0_dp*strength
               END DO
               resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
                                                2.0_dp*TARGET*strength*atom_coef(m)
            ELSE
               resp_env%matrix(atom_list_res(m), atom_list_res(m)) = &
                  resp_env%matrix(atom_list_res(m), atom_list_res(m)) + &
                  2.0_dp*strength
               resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
                                                2.0_dp*TARGET*strength
            END IF
         END DO
         DEALLOCATE (atom_list_res)
      END DO

      ! if heavies are restrained to zero, add these as well
      IF (resp_env%rheavies) THEN
         DO i = 1, natom
            CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, z=z)
            IF (z .NE. 1) THEN
               resp_env%matrix(i, i) = resp_env%matrix(i, i) + 2.0_dp*resp_env%rheavies_strength
            END IF
         END DO
      END IF

      !*** add the constraints
      ncons_v = 0
      ncons_v = ncons_v + natom

      ! REPEAT charges: treat the offset like a constraint
      IF (resp_env%use_repeat_method) THEN
         ncons_v = ncons_v + 1
         resp_env%matrix(1:natom, ncons_v) = resp_env%sum_vpot(1:natom)
         resp_env%matrix(ncons_v, 1:natom) = resp_env%sum_vpot(1:natom)
         resp_env%matrix(ncons_v, ncons_v) = 2.0_dp
         resp_env%rhs(ncons_v) = resp_env%sum_vhartree
      END IF

      ! total charge constraint
      IF (resp_env%itc) THEN
         ncons_v = ncons_v + 1
         resp_env%matrix(1:natom, ncons_v) = 1.0_dp
         resp_env%matrix(ncons_v, 1:natom) = 1.0_dp
         resp_env%rhs(ncons_v) = dft_control%charge
      END IF

      ! explicit constraints
      DO i = 1, resp_env%ncons_sec
         CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
         IF (.NOT. resp_env%equal_charges) THEN
            ncons_v = ncons_v + 1
            CALL section_vals_val_get(cons_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
            CALL section_vals_val_get(cons_section, "TARGET", i_rep_section=i, r_val=TARGET)
            CPASSERT(SIZE(atom_list_cons) == SIZE(atom_coef))
            DO m = 1, SIZE(atom_list_cons)
               resp_env%matrix(atom_list_cons(m), ncons_v) = atom_coef(m)
               resp_env%matrix(ncons_v, atom_list_cons(m)) = atom_coef(m)
            END DO
            resp_env%rhs(ncons_v) = TARGET
         ELSE
            my_atom_coef(1) = 1.0_dp
            my_atom_coef(2) = -1.0_dp
            DO k = 2, SIZE(atom_list_cons)
               ncons_v = ncons_v + 1
               resp_env%matrix(atom_list_cons(1), ncons_v) = my_atom_coef(1)
               resp_env%matrix(ncons_v, atom_list_cons(1)) = my_atom_coef(1)
               resp_env%matrix(atom_list_cons(k), ncons_v) = my_atom_coef(2)
               resp_env%matrix(ncons_v, atom_list_cons(k)) = my_atom_coef(2)
               resp_env%rhs(ncons_v) = 0.0_dp
            END DO
         END IF
         DEALLOCATE (atom_list_cons)
      END DO
      CALL timestop(handle)

   END SUBROUTINE add_restraints_and_constraints

! **************************************************************************************************
!> \brief print input information
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param rep_sys structure for repeating input sections defining fit points
!> \param my_per ...
! **************************************************************************************************
   SUBROUTINE print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(resp_p_type), DIMENSION(:), POINTER           :: rep_sys
      INTEGER, INTENT(IN)                                :: my_per

      CHARACTER(len=*), PARAMETER :: routineN = 'print_resp_parameter_info'

      CHARACTER(len=2)                                   :: symbol
      INTEGER                                            :: handle, i, ikind, kind_number, nkinds, &
                                                            output_unit
      REAL(KIND=dp)                                      :: conv, eta_conv
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_vals_type), POINTER                   :: input, resp_section

      CALL timeset(routineN, handle)
      NULLIFY (logger, input, resp_section)

      CALL get_qs_env(qs_env, &
                      input=input, &
                      atomic_kind_set=atomic_kind_set)
      resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
      logger => cp_get_default_logger()
      output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
                                         extension=".resp")
      nkinds = SIZE(atomic_kind_set)

      conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
      IF (.NOT. my_per == use_perd_none) THEN
         eta_conv = cp_unit_from_cp2k(resp_env%eta, "angstrom", power=-2)
      END IF

      IF (output_unit > 0) THEN
         WRITE (output_unit, '(/,1X,A,/)') "STARTING RESP FIT"
         IF (resp_env%use_repeat_method) THEN
            WRITE (output_unit, '(T3,A)') &
               "Fit the variance of the potential (REPEAT method)."
         END IF
         IF (.NOT. resp_env%equal_charges) THEN
            WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons_sec
         ELSE
            IF (resp_env%itc) THEN
               WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons - 1
            ELSE
               WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons
            END IF
         END IF
         WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit restraints: ", resp_env%nrest_sec
         WRITE (output_unit, '(T3,A,T80,A)') "Constrain total charge ", MERGE("T", "F", resp_env%itc)
         WRITE (output_unit, '(T3,A,T80,A)') "Restrain heavy atoms ", MERGE("T", "F", resp_env%rheavies)
         IF (resp_env%rheavies) THEN
            WRITE (output_unit, '(T3,A,T71,F10.6)') "Heavy atom restraint strength: ", &
               resp_env%rheavies_strength
         END IF
         WRITE (output_unit, '(T3,A,T66,3I5)') "Stride: ", resp_env%stride
         IF (resp_env%molecular_sys) THEN
            WRITE (output_unit, '(T3,A)') &
               "------------------------------------------------------------------------------"
            WRITE (output_unit, '(T3,A)') "Using sphere sampling"
            WRITE (output_unit, '(T3,A,T46,A,T66,A)') &
               "Element", "RMIN [angstrom]", "RMAX [angstrom]"
            DO ikind = 1, nkinds
               CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
                                    kind_number=kind_number, &
                                    element_symbol=symbol)
               WRITE (output_unit, '(T3,A,T51,F10.5,T71,F10.5)') &
                  symbol, &
                  resp_env%rmin_kind(kind_number)*conv, &
                  resp_env%rmax_kind(kind_number)*conv
            END DO
            IF (my_per == use_perd_none) THEN
               WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box min [angstrom]: ", resp_env%box_low(1:3)*conv
               WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box max [angstrom]: ", resp_env%box_hi(1:3)*conv
            END IF
            WRITE (output_unit, '(T3,A)') &
               "------------------------------------------------------------------------------"
         ELSE
            WRITE (output_unit, '(T3,A)') &
               "------------------------------------------------------------------------------"
            WRITE (output_unit, '(T3,A)') "Using slab sampling"
            WRITE (output_unit, '(2X,A,F10.5)') "Index of atoms defining the surface: "
            DO i = 1, SIZE(rep_sys)
               IF (i > 1 .AND. ALL(rep_sys(i)%p_resp%atom_surf_list == rep_sys(1)%p_resp%atom_surf_list)) EXIT
               WRITE (output_unit, '(7X,10I6)') rep_sys(i)%p_resp%atom_surf_list
            END DO
            DO i = 1, SIZE(rep_sys)
               IF (i > 1 .AND. ALL(rep_sys(i)%p_resp%range_surf == rep_sys(1)%p_resp%range_surf)) EXIT
               WRITE (output_unit, '(T3,A,T61,2F10.5)') &
                  "Range for sampling above the surface [angstrom]:", &
                  rep_sys(i)%p_resp%range_surf(1:2)*conv
            END DO
            DO i = 1, SIZE(rep_sys)
               IF (i > 1 .AND. rep_sys(i)%p_resp%length == rep_sys(1)%p_resp%length) EXIT
               WRITE (output_unit, '(T3,A,T71,F10.5)') "Length of sampling box above each"// &
                  " surface atom [angstrom]: ", rep_sys(i)%p_resp%length*conv
            END DO
            WRITE (output_unit, '(T3,A)') &
               "------------------------------------------------------------------------------"
         END IF
         IF (.NOT. my_per == use_perd_none) THEN
            WRITE (output_unit, '(T3,A,T71,F10.5)') "Width of Gaussian charge"// &
               " distribution [angstrom^-2]: ", eta_conv
         END IF
         CALL m_flush(output_unit)
      END IF
      CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
                                        "PRINT%PROGRAM_RUN_INFO")

      CALL timestop(handle)

   END SUBROUTINE print_resp_parameter_info

! **************************************************************************************************
!> \brief print RESP charges to an extra file or to the normal output file
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param output_runinfo ...
!> \param natom number of atoms
! **************************************************************************************************
   SUBROUTINE print_resp_charges(qs_env, resp_env, output_runinfo, natom)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      INTEGER, INTENT(IN)                                :: output_runinfo, natom

      CHARACTER(len=*), PARAMETER :: routineN = 'print_resp_charges'

      CHARACTER(LEN=default_path_length)                 :: filename
      INTEGER                                            :: handle, output_file
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(section_vals_type), POINTER                   :: input, print_key, resp_section

      CALL timeset(routineN, handle)

      NULLIFY (particle_set, qs_kind_set, input, logger, resp_section, print_key)

      CALL get_qs_env(qs_env, input=input, particle_set=particle_set, &
                      qs_kind_set=qs_kind_set)

      resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
      print_key => section_vals_get_subs_vals(resp_section, &
                                              "PRINT%RESP_CHARGES_TO_FILE")
      logger => cp_get_default_logger()

      IF (BTEST(cp_print_key_should_output(logger%iter_info, &
                                           resp_section, "PRINT%RESP_CHARGES_TO_FILE"), &
                cp_p_file)) THEN
         output_file = cp_print_key_unit_nr(logger, resp_section, &
                                            "PRINT%RESP_CHARGES_TO_FILE", &
                                            extension=".resp", &
                                            file_status="REPLACE", &
                                            file_action="WRITE", &
                                            file_form="FORMATTED")
         IF (output_file > 0) THEN
            filename = cp_print_key_generate_filename(logger, &
                                                      print_key, extension=".resp", &
                                                      my_local=.FALSE.)
            CALL print_atomic_charges(particle_set, qs_kind_set, output_file, title="RESP charges:", &
                                      atomic_charges=resp_env%rhs(1:natom))
            IF (output_runinfo > 0) WRITE (output_runinfo, '(2X,A,/)') "PRINTED RESP CHARGES TO FILE"
         END IF

         CALL cp_print_key_finished_output(output_file, logger, resp_section, &
                                           "PRINT%RESP_CHARGES_TO_FILE")
      ELSE
         CALL print_atomic_charges(particle_set, qs_kind_set, output_runinfo, title="RESP charges:", &
                                   atomic_charges=resp_env%rhs(1:natom))
      END IF

      CALL timestop(handle)

   END SUBROUTINE print_resp_charges

! **************************************************************************************************
!> \brief print potential generated by RESP charges to file
!> \param qs_env the qs environment
!> \param resp_env the resp environment
!> \param particles ...
!> \param natom number of atoms
!> \param output_runinfo ...
! **************************************************************************************************
   SUBROUTINE print_pot_from_resp_charges(qs_env, resp_env, particles, natom, output_runinfo)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(resp_type), POINTER                           :: resp_env
      TYPE(particle_list_type), POINTER                  :: particles
      INTEGER, INTENT(IN)                                :: natom, output_runinfo

      CHARACTER(len=*), PARAMETER :: routineN = 'print_pot_from_resp_charges'

      CHARACTER(LEN=default_path_length)                 :: my_pos_cube
      INTEGER                                            :: handle, ip, jx, jy, jz, unit_nr
      LOGICAL                                            :: append_cube, mpi_io
      REAL(KIND=dp)                                      :: dvol, normalize_factor, rms, rrms, &
                                                            sum_diff, sum_hartree, udvol
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(pw_c1d_gs_type)                               :: rho_resp, v_resp_gspace
      TYPE(pw_env_type), POINTER                         :: pw_env
      TYPE(pw_poisson_type), POINTER                     :: poisson_env
      TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      TYPE(pw_r3d_rs_type)                               :: aux_r, v_resp_rspace
      TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
      TYPE(section_vals_type), POINTER                   :: input, print_key, resp_section

      CALL timeset(routineN, handle)

      NULLIFY (auxbas_pw_pool, logger, pw_env, poisson_env, input, print_key, &
               para_env, resp_section, v_hartree_rspace)
      CALL get_qs_env(qs_env, &
                      input=input, &
                      para_env=para_env, &
                      pw_env=pw_env, &
                      v_hartree_rspace=v_hartree_rspace)
      normalize_factor = SQRT((resp_env%eta/pi)**3)
      resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
      print_key => section_vals_get_subs_vals(resp_section, &
                                              "PRINT%V_RESP_CUBE")
      logger => cp_get_default_logger()

      !*** calculate potential generated from RESP charges
      CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
                      poisson_env=poisson_env)

      CALL auxbas_pw_pool%create_pw(rho_resp)
      CALL auxbas_pw_pool%create_pw(v_resp_gspace)
      CALL auxbas_pw_pool%create_pw(v_resp_rspace)

      CALL pw_zero(rho_resp)
      CALL calculate_rho_resp_all(rho_resp, resp_env%rhs, natom, &
                                  resp_env%eta, qs_env)
      CALL pw_zero(v_resp_gspace)
      CALL pw_poisson_solve(poisson_env, rho_resp, &
                            vhartree=v_resp_gspace)
      CALL pw_zero(v_resp_rspace)
      CALL pw_transfer(v_resp_gspace, v_resp_rspace)
      dvol = v_resp_rspace%pw_grid%dvol
      CALL pw_scale(v_resp_rspace, dvol)
      CALL pw_scale(v_resp_rspace, -normalize_factor)
      ! REPEAT: correct for offset, take into account that potentials have reverse sign
      ! and are scaled by dvol
      IF (resp_env%use_repeat_method) THEN
         v_resp_rspace%array(:, :, :) = v_resp_rspace%array(:, :, :) - resp_env%offset*dvol
      END IF
      CALL v_resp_gspace%release()
      CALL rho_resp%release()

      !***now print the v_resp_rspace%pw to a cube file if requested
      IF (BTEST(cp_print_key_should_output(logger%iter_info, resp_section, &
                                           "PRINT%V_RESP_CUBE"), cp_p_file)) THEN
         CALL auxbas_pw_pool%create_pw(aux_r)
         append_cube = section_get_lval(resp_section, "PRINT%V_RESP_CUBE%APPEND")
         my_pos_cube = "REWIND"
         IF (append_cube) THEN
            my_pos_cube = "APPEND"
         END IF
         mpi_io = .TRUE.
         unit_nr = cp_print_key_unit_nr(logger, resp_section, &
                                        "PRINT%V_RESP_CUBE", &
                                        extension=".cube", &
                                        file_position=my_pos_cube, &
                                        mpi_io=mpi_io)
         udvol = 1.0_dp/dvol
         CALL pw_copy(v_resp_rspace, aux_r)
         CALL pw_scale(aux_r, udvol)
         CALL cp_pw_to_cube(aux_r, unit_nr, "RESP POTENTIAL", particles=particles, &
                            stride=section_get_ivals(resp_section, &
                                                     "PRINT%V_RESP_CUBE%STRIDE"), &
                            mpi_io=mpi_io)
         CALL cp_print_key_finished_output(unit_nr, logger, resp_section, &
                                           "PRINT%V_RESP_CUBE", mpi_io=mpi_io)
         CALL auxbas_pw_pool%give_back_pw(aux_r)
      END IF

      !*** RMS and RRMS
      sum_diff = 0.0_dp
      sum_hartree = 0.0_dp
      rms = 0.0_dp
      rrms = 0.0_dp
      DO ip = 1, resp_env%npoints_proc
         jx = resp_env%fitpoints(1, ip)
         jy = resp_env%fitpoints(2, ip)
         jz = resp_env%fitpoints(3, ip)
         sum_diff = sum_diff + (v_hartree_rspace%array(jx, jy, jz) - &
                                v_resp_rspace%array(jx, jy, jz))**2
         sum_hartree = sum_hartree + v_hartree_rspace%array(jx, jy, jz)**2
      END DO
      CALL para_env%sum(sum_diff)
      CALL para_env%sum(sum_hartree)
      rms = SQRT(sum_diff/resp_env%npoints)
      rrms = SQRT(sum_diff/sum_hartree)
      IF (output_runinfo > 0) THEN
         WRITE (output_runinfo, '(2X,A,T69,ES12.5)') "Root-mean-square (RMS) "// &
            "error of RESP fit:", rms
         WRITE (output_runinfo, '(2X,A,T69,ES12.5,/)') "Relative root-mean-square "// &
            "(RRMS) error of RESP fit:", rrms
      END IF

      CALL v_resp_rspace%release()

      CALL timestop(handle)

   END SUBROUTINE print_pot_from_resp_charges

END MODULE qs_resp
